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Nanoflows through disordered media: a joint Lattice Boltzmann 
and Molecular Dynamics investigation 
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Abstract. - We investigate nanoflows through dilute disordered media by means of joint lattice 
Boltzmann (LB) and molecular dynamics (MD) simulations — when the size of the obstacles is 
comparable to the size of the flowing particles — for randomly located spheres and for a correlated 
particle-gel. In both cases at sufficiently low solid fraction, $ < 0.01, LB and MD provide similar 
values of the permeability. However, for $ > 0.01, MD shows that molecular size effects lead 
to a decrease of the permeability, as compared to the Navier-Stokes predictions. For gels, the 
simulations highlights a surplus of permeability, which can be accommodated within a rescaling 
of the effective radius of the gel monomers. 



Introduction. — Flow phenomena in disordered me- 
dia are a subject of great theoretical and practical interest 
[HE] • in case 01 a nm d streaming through a random, low- 
density porous matrix, descriptions based on continuum 
hydrodynamics have been provided [3 8 . Continuum the- 
ories are expected to hold on macroscopic scales and thus 
might be non-applicable to fluid flow through microfluidic 
devices or through micro-gel matrices where the typical 
size of the pores is at nanometer length scales and below. 

However, as shown in several simulation studies (see 
e.g. Refs. [51110]). hydrodynamics often holds down to the 
molecular scale, as far as simple steady state flows of dense 
liquids are concerned. Similar conclusions have also been 
reached recently for the non-trivial case of microflows over 
super- hydrophobic surfaces QT]. However, in particular 
for the case of nanoflow in disordered media, the question 
of whether/to what extent continuum theory is applica- 
ble to flows at the nanoscale in porous materials, remains 
open to this day. Indeed, previous studies on flow phenom- 
ena in disordered media have employed only mesoscopic 
or macroscopic simulation methods, such as finite clement 
schemes [EJUS], the LB method (T4l!21j . and smoothed 
particle dynamics [22 ] f23 ] . In this work, we address such a 
question using a combination of Lattice-Boltzmann (LB) 
and molecular dynamics (MD) simulations. While LB is 
used as an effective Navier-Stokes equation solver, the MD 



simulations are performed to solve Newton's equations of 
motion for a three-dimensional system of soft spheres. As 
porous media, we consider non-overlapping random ar- 
rangements of particles, as well as particle gel networks at 
different packing fractions. In both cases, flows through 
bulk porous media and through porous media confined be- 
tween parallel plates are considered. First, we show that 
a quantitative mapping between LB and MD can be es- 
tablished. Second, we study to what extent continuum 
theory is applicable at the molecular scale. This is par- 
ticularly interesting for gel networks, since they introduce 
long-ranged structural correlations which are not present 
in the random arrangement of obstacles. For this case, 
we quantify deviations to predictions of the theory. In- 
terestingly, the theory is rcnormalizable, i.e. the observed 
deviations due to an excess of permeability can be reab- 
sorbed within a readjustment of the effective radius of the 
gel monomers. 

Theory. — For low Reynolds number flow, Darcy 
[24] first established empirically a linear relation between 
the average volumetric flow velocity through unit cross- 
sectional area, ({/), and the pressure gradient Vp of the 
fluid across the porous medium, (v) = — (k/rj)Wp (with 
r\ the shear viscosity of the fluid and k the permeability). 
Note that the measurement or calculation of (v) implies an 
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average over different reaiizations of the porous medium. 
In case of high dilution of spheres of radius R as a porous 
medium, Darcy's law reduces essentially to Stokes' law, 
with a permeability k = 2R 2 /%4>, [Hfll2"5] 

Brinkman's theory is based on the stationary Navier- 
Stokes equations r)V 2 v — Vp + F = and V ■ V = 0, 
with v the flow velocity and F an external force per unit 
volume acting on the fluid. By setting v = (v) and by 
expressing the external force F through Darcy's law, F = 
—rjv/k, Brinkman's equation of motion [3] is obtained: 
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This equation describes the porous matrix as an effective 
medium that exerts a friction on the fluid. The substitu- 
tion of the external force F by Darcy's law is expected to 
be valid only if the packing fraction of the porous medium 
is sufficiently small (sec below). 

From Eq. ([1]), one can derive an explicit formula for 
the flow between two parallel plates in presence of a 
porous medium. Consider a gravitational force field pg 
in x-direction and two parallel plates at z = —L/2 and 
z = L/2. The porous medium between the plates is 
represented by a random matrix of fixed non-overlapping 
spheres. Then, the velocity profile is given by: 
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with v = r\j p the kinematic velocity of the fluid. 

For a dilute collection of non-overlapping spheres, dif- 
ferent expressions for the permeabiliy as a function of the 
volume fraction, $ = 47ri? 3 n/3, of the porous medium 
have been investigated [3HZ]- A simple expression for the 
permeability over the entire range of volume fractions has 
been proposed by van der Hoef et al. [7] by fitting both 
LB [SJH] and multipole expansion data [5] 
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Note that the square-root of the permeability, y/k, de- 
scribes the screening length of the flow field due to the 
interaction with the porous medium. 

Equations ([2j and ([3]) are the main results which our 
simulations on nanoscopic scales will be compared against. 

Methods and simulation details. We use a stan- 
dard Lattice Boltzmann model with single-step relaxation 
term [26H28] . No-slip boundary conditions at solid sur- 
faces are implemented via standard bounce-back collision 
rules [29]. Fully developed periodic flows are generated 
by pressure boundary conditions [30) . The size of the LB 
D3Q19 lattice is (192 x 48 x 48) in lattice units. When 
present, slit walls are placed parallel to the xy plane with 
surfaces at locations z = 4 and z = 43. 




Fig. 1: Two snapshots at $ = 0.1 with random medium (a) 
and gel (b). Porous media are coded in red, the binary mixture 
fluid is coded in yellow and green. The size of the fluid particles 
is not at scale to improve visibility. 



For the MD fluid, we choose a binary model system as 
proposed by Hedges et al. [31] for which, also at low tem- 
peratures, crystallization is not a problem and thus this 
model can be used in forthcoming studies on glassforming 
fluids in porous media. 

The MD fluid is composed of a 50:50 mixture of A 
and B type particles interacting with a WCA poten- 
tial, V a p(r) = {ie a i3[{r/(T al3 y 12 - {r/a al3 y 6 } + e afi } for 
r < r_ = 2 1 / 6 er a/ 3 (zero else). The parameters are 
£AA = £AB = e B B = 1 and (TAA = 1.0, (TAB = 11/12, 
cbb = 5/6 [21]. In the following, length and energy scales 
are measured in units of a = oaa and e = caa, respec- 
tively; temperature is in units of the potential depth e. 
The masses of both A and B particles are taken equal, 
m = niA = tub = 1. The equations of motion are in- 
tegrated with the velocity form of the Verlet algorithm 
using a time step 5t = 10~ 3 in units of to = [ma 2 /e] 1 ^ 2 . 
Thcrmalization at the constant temperature T = 5.0 is 
obtained with a Lowe thermostat |32j . which provides lo- 
cal momentum conservation and thus it preserves the cor- 
rect hydrodynamic behavior of the fluid. The size of the 
MD simulation box is (32 x 8 x 8) in units of a. The 
total density of all systems is set to p = 1.3. First the 
porous material is introduced in the simulation box which 
is then filled with the appropriate number of fluid particles 
to reach the desired density. A pressure drop is applied 
along the x direction by adding a gravitational field (g) 
on the fluid particles, whose intensity (g in units of ct/Vq) 
is chosen such that a linear response of the system is pro- 
vided. The kinematic viscosity v of the fluid has been 
calculated aside from separate Poiseuille-flow simulations, 
which yields v = 24 a 2 /to. 

To simulate slit walls (whenever present), the system is 
first equilibrated in a (32 x 8 x 10.2) a box and then all 
fluid particles within a distance z w = 1.2a from the planes 
z = and z = 10.2cr, are labelled as wall particles. Wall 
particles retain the same interactions as fluid particles, 
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but their position is kept fixed. In this way we account for 
the walls roughness and avoid layering effects. To avoid 
diffusion of fluid particles in the wall region, an external 
potential (previously defined as Vaa) is added along two 
planes parallel to the xy plane, at a distance z p — 2 1 / 6 er 
from the slit walls. In this way, fluid particles feel the 
external potential only when they start diffusing into the 
walls. 

To compare LB and MD, units have to be scaled ap- 
propriately. One of our goals is to assess the validity of 
the LB predictions for microscopic flows which can be re- 
solved through MD. As shown recently [10], in order to 
recover quantitative agreement with MD results, the LB 
simulation must be taken down to microscopic resolution, 
i.e. fractions of the range of molecular interactions. The 
space conversion proceeds as follows. In LB simulations 
the lattice spacing, Ax is set to unity, while in MD simula- 
tions the unit of length is fixed by the parameter a. To re- 
solve fractions of the interaction potential we set Ax < a, 
specifically we choose a = 6Ax. This fixes the conversion 
of space units. It also fixes the radius of the LB spheres 
as R = a/2 = 3 Ax. In the rest of the letter we adopt 
c = &aa = 1 as our unit of length. 

The time conversion is determined from kinematic vis- 
cosity v. The kinematic viscosity in the LB simulation is 

A 2 

given by z^lb = z^lb -£p , with z-'lb = 1/6, while for MD we 

~ a 2 ~ 

have fMD = ^md Atltu '■ w ^ n ^md given by the Poisscuillc 
flow comparison. By imposing z^lb = ^md and remember- 
ing the space conversion factor we obtain for time scales 



The Reynolds number is Re =< v > R/v, where < 
v > is the spatially averaged velocity, R is the radius of 
the obstacles and v is the kinematic viscosity. We use 
very low external fields (g — 0.05 in units of a/At 2 MD ), so 
that Re < 10 -3 and simulations can be considered under 
effectively zero-Reynolds number conditions. 

MD and LB simulations are best compared in terms of 
dimensionless quantities, such as the normalized perme- 
ability fc/fco , where kg is the permeability of a single sphere 
(as given by Stokes law). Other dimensionless quantities 
used in the present work are the dimensionless position 
z/L, with L the width of the slit, and the dimensionless 
velocity profile u = vv / (gR 2 ). 

Simulations of fluid flow through two types of obstacles 
are considered: random media and gel media. 

Random media are modeled as a collection of non- 
overlapping spheres of radius R. We average over 50 in- 
dependent random configuration for each volume fraction 
$ considered, in the range $ = 0.008 — 0.27. 

Gel networks are characteristic random structures with 
long range spatial correlations. A gel is usually made of 
a network of polymer strands or from self-assembled col- 
loidal particles. The number of bonds is so high that it is 
always possible to move from one gel-forming monomer to 
another without ever leaving the network. This interlinked 
structure confers the gel its peculiar properties, sharing 



characteristics of both liquids and solids. It behaves like 
a liquid since it can be made up primarily of fluid and 
allow both diffusive and convective transport through its 
volume. On the other hand, a gel is also able to support a 
shear stress and behave elastically, acting like a solid. In 
the present study, we neglect the elastic behaviour of gels, 
keeping the monomers fixed and taking advantage of the 
rigid framework through which mass transport can occur. 
Gel structures are obtained through equilibrium MD sim- 
ulations of Patchy particles. Patchy particles are a class of 
short-ranged valence-limited particles which can reach low 
temperatures without encountering the gas-liquid phase 
separation [33) . The corresponding low T arrested states 
are in fact equilibrium gels which we use for the present 
study. We follow the procedure described in [53], for net- 
works with average valence 2.25, equilibrated until all par- 
ticles belong to the same spanning cluster. Here, we gen- 
crated 50 independent gel configurations at each of the 
packing fractions $ = 0.03, 0.05, 0.075, 0.1, 0.15, and 0.2. 
Figure [T] shows snapshots of MD configurations with the 
random medium (a) and the gel (b) at $ = 0.1. 

Results. — All obstacles considered in the present pa- 
per are collections of spheres whose hydrodynamic radius 
is estimated by measuring the drag force on a single sphere 
with periodic boundary conditions and comparing the re- 
sult with the theoretical expression given by Hashimoto 
[35] . For LB, this procedure leads to the an hydrodynamic 
radius of R = 3.05, slighty larger than the nominal radius 
(in agreement with previous studies [2H])- For MD we 
obtain a value of R = 0.44, assuming slip boundary condi- 
tions. The same value, R = 0.44 can be estimated from the 
interaction potential for the MD particles: at 2 • R = 0.88 
the potential energy between two particles, separated by a 
distance 2 • R = 0.88, is approximately equal to ksT. We 
conclude that for our atomistic fluid slip boundary condi- 
tions provide a consistent description of the interactions 
between obstacles and fluid particles. We remind that the 
difference between stick and slip boundary conditions is 
included in the Stokes expression: fco = 2i? 2 /(9$) (stick) 
and fc = i? 2 /(3$) (slip). 

Figure [2] shows the results for the flow through random 
porous media (without slit walls) for both MD and LB 
simulations. LB results confirm Eq. ([3]), except for the 
lowest volume fractions, where the low-density result of 
result of Kim and Russel [4] holds. Note that the agree- 
ment with Brinkman theory is also good for low volume 
fractions, say $ < 0.08. The grid-independence of LB re- 
sults was checked by doubling the resolution [a = 12 Ax). 
MD data shows instead significant deviations from the hy- 
drodynamic solution. While for $ < 0.01 the agreement 
between MD and LB seems to hold (within statistical ac- 
curacy) , at higher volume fractions the MD results show 
visible under-deviations. These can be interpreted as gen- 
uine atomistic effects, where the reduced permeability is 
due to the finite size of the molecules. Indeed, at $ = 0.01, 
the ratio d/a between the average intermolecular distance 
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Fig. 2: Dimensionless permeability as a function of volume 
fraction $ for the random sphere matrix from MD simulations 
(circles) and LB simulations (squares). The LB data is well 
described by the theoretical prediction of Eq. (J3j (continuous 
line), except at low volume fractions where the agreement holds 
with the asymptotically exact result of Kim and Russel [4] 
(dashed line). For comparison also the result of the Brinkman 
theory [3j (dotted line) is shown. MD results consistently de- 
viate from hydrodynamic predictions for <J> > 0.01. 

and the molecular effective diameter is of the order of 
d/cr ru 4, whereas for $ = 0.1, one has dja ~ 2. 

We have also investigated the effect of slit walls on the 
permeability k. To this end, we have determined the veloc- 
ity profiles v x (z), averaged over all crossflow yz sections. 
Inspection of the LB and MD velocity profiles in Fig. © 
for random media at different values of $ shows excel- 
lent agreement with the prediction of Eq. ((SJ for LB. In 
contrast to that, the MD profiles show a different pattern 
from the hydrodynamic solution: again, this indicates the 
differences between hydrodynamic and atomistic flow for 
$ > 0.01. 

We can conclude that both global (permeability) and lo- 
cal (velocity profiles) properties of an atomistic flow show 
consistent deviations from the hydrodynamic solution also 
at relatively low volume fractions. The results obtained 
show the limits of applicability of continuum approaches 
when dealing with microscopic and structured fluids. 

We now turn to the study of the effects of the correlation 
between the obstacles' positions on the flow properties for 
both the hydrodynamic solution (LB) and the microscopic 
flow (MD). Differently from random spheres, gel media is 
characterized by short range correlations between the par- 
ticles' positions. Figure @] reports the results for the gel 
media for both MD and LB simulations. This figure in- 
dicates that the permeability of the gel is always higher 
than that of the random media for both hydrodynamic 
and microscopic flows. Also for the gel, the MD perme- 
ability always underestimates the LB solution at packing 
fractions $ > 0.01. 

Next, we show that the increased permeability induced 
by the correlations between the obstacle particles can be 
taken into account by simple renormalization of the theory 
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Fig. 3: Comparison of LB (squares) and MD (circles) for the 
random medium at different volume fractions. Full lines are 
the theoretical profiles, Eq. @. Agreement between LB data 
and theoretical predictions is perfect for all volume fractions 
investigated. MD solutions display instead a different pattern. 



for random obstacles. We first observe that the permeabil- 
ity can be formally written as k = i? 2 /($), where /($) is 
a dimensionless function of 4>. Given the proper expres- 
sion for the random sphere case [Eq. (J3| for the LB data] , 
a formal generalization for the gel case can be obtained by 
introducing an effective parameter a($), 

k = a($)i? 2 /($) ■ (4) 

The introduction of the parameter a($) formally cor- 
responds to the definition of an effective Stokes radius, 
i?cff($) = y/aR, which depends on density. The inset of 
Fig. [2] shows the value of the parameter a calculated both 
for LB (squares) and MD (circles) simulations. For both 
types of simulations the effective Stokes radius remains al- 
most constant at all considered volume fractions, so that 
we can conclude that the gel structure results in an in- 
crease of the bare permeability by about 50%, which can 
be also interpreted as an increase of the effective Stokes 
radius by a factor 1.25. 

Summary and conclusions. — We have investigated 
nanoflows through disordered media by means of joint LB 
and MD simulations. For random media at $ < 0.01, 
both LB and MD provide similar values of the perme- 
ability, confirming that the hydrodynamic approach holds 
down to the nanoscale at sufficiently low volume fractions. 
For higher volume fractions, the MD simulations reveal 
undcr-departures from the hydrodynamic solution. Since 
LB simulations still agree with the Eq. [3] prediction, these 
under-departures are most naturally interpreted as gen- 
uine atomistic effects, i.e. breakdown of the continuum 
hydrodynamic hypothesis at the nanoscale. The same 
conclusions result from the study of the velocity profiles 
within solid slabs filled with random media, showing that 
the MD flow is qualitatively different from LB predictions. 
We have then explored the effects of the correlations in the 
obstacles positions by studying the flow through a corre- 
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Fig. 4: Dimensionless permeability as a function of volume 
fraction $ for the gel matrix from MD simulations (circles) 
and LB simulations (squares). The dashed line represents the 
results for the random sphere matrix from Fig. [2] The inset 
shows the a parameter defined in Eq. [4] for MD (circles) and 
LB (squares). 



lated medium, i.e. gels. We have shown that gels exhibit 
a surplus of permeability. This surplus, about 50%, can 
be reinterpreted as an increase of the effective radius of 
the gel monomers, thereby indicating that, at least for the 
cases explored in this work, the hydrodynamic solution 
appears to be renormalizable. 
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